**Religiosity Chapter**

set more off
clear
capture log close
cd "~/Documents/Projects/Religiosity/ Chapter"
log using relig_chapter.log, replace

local data_dir = "~/Documents/Projects/Data"

// WVS/EVS (Five waves)
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
use "`data_dir'/WVS/xwvs+.dta", clear
keep s003 s025 country year f028  
append using "`data_dir'/WVS/Waves 1-5/new_wvs1-5.dta", nolabel  ///
	keep(s003 s025 country year f028)


*church
gen church8=9-f028
label var church8 "Church attendance, 8 categories"

gen church2=church8
recode church2 1/6=0 7/8=1
label var church2 "Attends church weekly or more often"

encode country, gen(ccode) label(ccode)
collapse s003 year ccode church2, by(s025)
label val ccode ccode
decode ccode, gen(country)

replace church2=church2*100

replace country="United Kingdom" if country=="Great Britain"

drop s025 s003 ccode
order country year 
sort country year
save wvs_evs.dta, replace

// EVS 2008
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
use "`data_dir'/EVS/2008/ZA4800_F1.dta", clear 
keep country cntry_y v109

*church
gen churchyr=v109
recode churchyr 1=104 2=52 3=12 4=2 5=1 6=.5 7=0
label var churchyr "Annual church attendance"

gen church2=churchyr
replace church2 = (church2 >= 52) if !mi(church2)
label var church2 "Attends church weekly or more often"

decode country, gen(ccode0) 
encode ccode0, gen(ccode) label(ccode)
collapse ccode church*, by(cntry_y)
label val ccode ccode
decode ccode, gen(country)

replace church2=church2*100

replace country="United Kingdom" if country=="Great Britain"

drop cntry_y ccode
gen year = 2008
order country year 
sort country year
save evs2008.dta, replace


// ISSP Religion (1991, 1998, 2008)
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
use "`data_dir'/ISSP/Religion/Religion I/ZA2150_F1.dta", clear


*church
gen churchyr = v107 if !mi(v107)
recode churchyr 1=52 2=30 3=12 4=3 5=1 6=0
label var churchyr "Annual church attendance"

gen church2=churchyr
replace church2 = (church2 >= 52) if !mi(church2)
label var church2 "Attends church weekly or more often"

recode v3 2=1
collapse church*, by(v3)

replace church2 = church2*100

gen country=""
replace country="Germany" if v3==1
replace country="United Kingdom" if v3==3
replace country="United States" if v3==5
replace country="Hungary" if v3==6
replace country="Netherlands" if v3==7
replace country="Italy" if v3==8
replace country="Ireland" if v3==9
replace country="Norway" if v3==10
replace country="Austria" if v3==11
replace country="Slovenia" if v3==12
replace country="Poland" if v3==13
replace country="Israel" if v3==14
replace country="Philippines" if v3==15
replace country="New Zealand" if v3==16
replace country="Australia" if v3==18
drop if country==""
drop v3

gen year=1991
replace year=1990 if country=="Italy"
replace year=1993 if country=="Austria"

order country year
sort country year
save issp_rel1.dta, replace

//1998
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
use "`data_dir'/ISSP/Religion/Religion II/ZA3190_F1.dta", clear

*church
gen churchyr = v218 if !mi(v218)
recode churchyr 1=52 2=30 3=12 4=3 5=1 6=0
label var churchyr "Annual church attendance"

gen church2=churchyr
replace church2 = (church2 >= 52) if !mi(church2)
label var church2 "Attends church weekly or more often"

recode v3 3=2
collapse church*, by(v3)
decode v3, gen(country)
replace country=regexs(1) if regexm(country, ".*-(.*)")
replace country="Israel" if country=="Israel Jews Arabs"
replace country="United Kingdom" if country=="Great Britain"
replace country="Germany" if country=="West"
drop v3

replace church2 = church2*100

gen year=1998
replace year=1999 if country=="Czech Republic" | country=="Austria" | country=="Poland"  ///
	| country=="Portugal" | country=="Switzerland"
replace year=2000 if country=="Canada"

order country year
sort country year
save issp_rel2.dta, replace

//2008
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
use "`data_dir'/ISSP/Religion/Religion III/ZA4950_F1.dta", clear

*church
gen churchyr = ATTD_E if !mi(ATTD_E)
recode churchyr 1=0 2=.5 3=1.5 4=3 5=12 6=30 7=45 8=52 9=104
label var churchyr "Annual church attendance"

gen church2=churchyr
replace church2 = (church2 >= 52) if !mi(church2)
label var church2 "Attends church weekly or more often"

collapse church*, by(V5)
decode V5, gen(country)
replace country=regexs(1) if regexm(country, ".*-(.*)")
replace country="Korea, Republic of" if country=="South Korea"
replace country="United Kingdom" if country=="Great Britain and/or United Kingdom"
replace country="Slovak Republic" if country=="Slovakia"
drop V5

replace church2 = church2*100

gen year=2008
replace year=2009 if country=="Australia" | country=="Croatia" | country=="Denmark" | country=="Israel"  ///
	| country=="Italy" | country=="Latvia" | country=="Mexico" | country=="Latvia"  ///
	| country=="Portugal" | country=="Russia" | country=="Slovenia" | country=="Taiwan"
replace year=2010 if country=="Poland"

order country year
sort country year
save issp_rel3.dta, replace



// ISSP Family (1994)
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
use "`data_dir'/ISSP/Family/Family II/ZA2620_F1.dta", clear

*church
gen churchyr = v220 if !mi(v220)
recode churchyr 1=52 2=30 3=12 4=3 5=1 6=0
label var churchyr "Annual church attendance"

gen church2=churchyr
replace church2 = (church2 >= 52) if !mi(church2)
label var church2 "Attends church weekly or more often"

recode v3 3=2
collapse church*, by(v3)
decode v3, gen(country)
replace country=regexs(1) if regexm(country, ".*- ?(.*)")
replace country="United Kingdom" if country=="Great Britain"
replace country="Germany" if country=="West"
drop v3

replace church2 = church2*100

gen year=1994
replace year=1995 if country=="Bulgaria" | country=="Austria"
replace year=1993 if country=="Canada" | country=="Slovenia"

order country year
sort country year
save issp_fam2.dta, replace


// Roper
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
insheet using "`data_dir'/Roper Abortion/roper.csv", comma clear
keep country year church2
keep if church2~=.
sort country year
save roper.dta, replace


// Other surveys without microdata
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
insheet using "`data_dir'/Roper Abortion/other.csv", comma clear
keep country year church2
keep if church2~=.
sort country year
save other.dta, replace


// Germany ALLBUS GSS
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
use "`data_dir'/GSS/Germany/ZA4572_GCUM.dta", clear

gen year = v2

*church
gen churchyr = v544
recode churchyr 1=104 2=52 3=24 4=3 5=1 6=0 7=. 9=. 0=.
label var churchyr "Annual church attendance"

gen church2=churchyr
replace church2 = (church2 >= 52) if !mi(church2)
label var church2 "Attends church weekly or more often"

collapse church2, by(year)

replace church2=church2*100

gen country="Germany"
order country year
sort country year
save de_gss.dta, replace


// Ireland Elections (2002-2007)
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
use "`data_dir'/Country Surveys/Ireland/INESLong_Beta.dta", clear

gen churchyr = v0964 if v0964<9
recode churchyr 1=104 2=52 3=30 4=12 6=3 7=1 8=.5 9=0

gen church2 = (churchyr>=52) if churchyr~=.

gen year = ines
drop if year == .
collapse church2, by(year)
replace church2 = church2*100
gen country="Ireland"
order country year
save ireland.dta, replace


// Poland Elections (1997, 2000, 2001)
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
use "`data_dir'/Country Surveys/Poland/1997/ZA4333_F1.dta", clear

gen churchyr = churchat if churchat<9
recode churchyr 1=0 2=.5 3=1 4=3 5=12 6=30 7=52 8=104

gen church2 = (churchyr>=52) if churchyr~=.

gen year = 1997
collapse church2, by(year)
replace church2 = church2*100

gen country = "Poland"
order country year
save poland97.dta, replace

//2001
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
use "`data_dir'/Country Surveys/Poland/2001/ZA4335_F1.dta", clear

gen churchyr = chrat365 if chrat365<9
recode churchyr 1=0 2=.5 3=1 4=3 5=12 6=30 7=52 8=104

gen church2 = (churchyr>=52) if churchyr~=.

gen year = 2001
collapse church*, by(year)
replace church2 = church2*100

gen country = "Poland"
order country year
save poland01.dta, replace

use poland97.dta, clear
append using poland01.dta
save poland.dta, replace


// U.S. American National Election Studies, 1948-2004
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
use "`data_dir'/ANES/1948-2004/08475-0001-Data.dta", clear

gen year=VCF0004

gen church2=(VCF0130==1) if year>1970

collapse church2, by(year)

egen c = rowmean(church2)
drop if c==.
drop c

local new = _N + 1
set obs `new'
recode year .=2006
recode church2 . = 0.248 if year==2006	/*ANES 2006 Pilot Study*/

local new = _N + 1
set obs `new'
recode year .=2008
recode church2 . = 0.244 if year==2008	/*ANES 2008 Pre- and Post-*/

gen country = "United States"

replace church2 = church2*100

save us_anes.dta, replace


// U.S. GSS, 1972-2008
if "`data_dir'" == "" local data_dir = "~/Documents/Projects/Data"
use "`data_dir'/GSS/U.S./25962-0001-Data.dta", clear

foreach v of varlist _all {
	rename `v' `=lower("`v'")'
	*"
}

*church
gen churchyr=attend
recode churchyr 0=0 1=.5 2=1.5 3=3 4=12 5=30 6=44 7=52 8=104
label var churchyr "Annual church attendance"

gen church2=churchyr
replace church2 = (church2 >= 52) if !mi(church2)
label var church2 "Attends church weekly or more often"


collapse church2, by(year)

replace church2=church2*100

gen country="United States"
order country year
sort country year
save us_gss.dta, replace

// Eurbarometer Trend
use "~/Documents/Projects/Data/EB Trend/EBtrend_2.dta", clear
gen church2 = (church>=52) & church~=.
rename church churchyr
collapse church2 churchyr, by(country year)
replace church2=church2*100
drop if church2==0
drop if year>=1995
save eb_trend.dta, replace


// European Social Survey Cumulative
use "~/Documents/Projects/Data/ESS/ESS1-4/ess1-4.dta", clear
gen church2 = (rlgatnd<4) & rlgatnd~=.
gen churchyr = rlgatnd
recode churchyr 1=365 2=104 3=52 4=24 5=2 6=.5 7=0
collapse church2 churchyr, by(country year)
replace church2=church2*100
save ess_cum.dta, replace


// International Social Mobility and Politics
use "~/Documents/Projects/Data/ISMP/p1145/p1145.dta", clear
gen church2 = (religatt==1) if religatt>0
kountry country, from(iso3c)
drop country
rename NAMES_STD country
replace country="Denmark" if country=="den"
replace country="United Kingdom" if country=="eng"
replace country="Ireland" if country=="ire"
replace country="Netherlands" if country=="net"
replace country="Switzerland" if country=="swi"
collapse church2, by(country year)
drop if church2==.
replace church2=church2*100
decode year, gen(yr)
drop year
destring yr, gen(year)
drop yr
keep if year>1962
save ismp.dta, replace


// ISSP
*Role of Government 1: 1985
use "~/Documents/Projects/Data/ISSP/Role of Gov't/Role of Gov't I/issp85_0.dta"
gen churchyr = v133
recode churchyr 1=52 2=24 3=3 4=1 5=0
gen church2 = (churchyr>=52) & !mi(churchyr)
collapse church2 churchyr, by(country year)
replace church2=church2*100
drop if mi(churchyr)
save issp_gov1.dta, replace


*Role of Government 2: 1990
use "~/Documents/Projects/Data/ISSP/Role of Gov't/Role of Gov't II/issp90_0.dta", clear
drop if v3==5
replace country="Germany" if v3==3
gen churchyr = v89
recode churchyr 1=52 2=30 3=12 4=3 5=1 6=0
gen church2 = (churchyr>=52) & !mi(churchyr)
collapse church2 churchyr, by(country year)
replace church2=church2*100
drop if mi(churchyr)
save issp_gov2.dta, replace


*National Identity 1: 1995
use "~/Documents/Projects/Data/ISSP/National Identity/National Identity I/issp95_0.dta", clear
gen churchyr = v266
recode churchyr 1=52 2=30 3=12 4=3 5=1 6=0
gen church2 = (churchyr>=52) & !mi(churchyr)
collapse church2 churchyr, by(country year)
replace church2=church2*100
drop if mi(churchyr)
save issp_nat1.dta, replace


* Role of Government 3: 1996
use "~/Documents/Projects/Data/ISSP/Role of Gov't/Role of Gov't III/issp96_0.dta", clear
gen churchyr = v220
recode churchyr 1=52 2=30 3=12 4=3 5=1 6=0
gen church2 = (churchyr>=52) & !mi(churchyr)
collapse church2 churchyr, by(country year)
replace church2=church2*100
drop if mi(churchyr)
save issp_gov3.dta, replace


*Social Networks 2: 2001
use "~/Documents/Projects/Data/ISSP/Social Networks/Social Networks II/issp01_0.dta", clear 
gen churchyr = attend if attend<9
recode churchyr 1=52 2=30 3=12 4=3 5=1 6=0 7=1
gen church2 = (churchyr>=52) & !mi(churchyr)
collapse church2 churchyr, by(country year)
replace church2=church2*100
drop if mi(churchyr)
save issp_net2.dta, replace


*National Identity 2: 2003
use "~/Documents/Projects/Data/ISSP/National Identity/National Identity II/issp03_0.dta", clear
gen churchyr = attend if attend<9
recode churchyr 1=104 2=52 3=30 4=12 5=3 6=1 7=.5 8=0
gen church2 = (churchyr>=52) & !mi(churchyr)
collapse church2 churchyr, by(country year)
replace church2=church2*100
drop if mi(churchyr)
save issp_nat2.dta, replace


*Citizenship 1: 2004
use "~/Documents/Projects/Data/ISSP/Citizenship/Citizenship I/issp04_0.dta", clear
gen churchyr = v300 if v300<9
recode churchyr 1=104 2=52 3=30 4=12 5=3 6=1 7=.5 8=0
gen church2 = (churchyr>=52) & !mi(churchyr)
collapse church2 churchyr, by(country year)
replace church2=church2*100
drop if mi(churchyr)
save issp_cit1.dta, replace


*Role of Government 4: 2006
use "~/Documents/Projects/Data/ISSP/Role of Gov't/Role of Gov't IV/issp06_0.dta", clear
gen churchyr = ATTEND if ATTEND<9
recode churchyr 1=104 2=52 3=30 4=12 5=3 6=1 7=.5 8=0
gen church2 = (churchyr>=52) & !mi(churchyr)
collapse church2 churchyr, by(country year)
replace church2=church2*100
drop if mi(churchyr)
save issp_gov4.dta, replace


// Assemble data
use wvs_evs.dta, clear
append using evs2008.dta
append using issp_rel1.dta
append using issp_rel2.dta
append using issp_rel3.dta
append using issp_fam2.dta
append using roper.dta
append using other.dta
append using de_gss.dta
append using ireland.dta
append using poland.dta
append using us_anes.dta
append using us_gss.dta
append using eb_trend.dta
append using ess_cum.dta
append using ismp.dta
append using issp_gov1.dta
append using issp_gov2.dta
append using issp_nat1.dta
append using issp_gov3.dta
append using issp_net2.dta
append using issp_nat2.dta
append using issp_cit1.dta
append using issp_gov4.dta

keep country year church2
drop if mi(year)
drop if country=="Andorra" | country=="Northern Ireland"
drop if country=="Australia" & year==1979
drop if country=="South Africa" & year==2008
replace country="Bosnia and Herzegovina" if country=="Bosnia Herzegovina"
replace country="Russian Federation" if country=="Russia"
replace country="Macedonia, FYR" if country=="Macedonia"
replace country="United Kingdom" if country=="Great Britain"
replace country="Slovak Republic" if country=="Slovakia"
replace country="Korea, Republic of" if country=="South Korea"
replace country="Yugoslavia" if country=="Serbia and Montenegro" & year<2001
drop if church2==.

encode country, gen(ccode) label(ccode)
collapse church2, by(ccode year)
label val ccode ccode
decode ccode, gen(country)
drop ccode
order country year

sort country year
save trend00.dta, replace

use trend00.dta, clear
bysort country: gen lag=year[_n+1]-year
drop if lag>5 & !mi(lag)
drop lag
bysort country: gen lag=year[_n+1]-year
drop if lag>5 & !mi(lag)
drop lag
bysort country: egen ok1= count(church2)
drop if ok1<7
drop ok1

encode country, gen(cc) label(cc)
tsset cc year
tsfill
sort cc year
drop country
decode cc, gen(country)
order country year
save trend0.dta, replace

use trend0.dta, clear
tsset cc year
order country year

sort country year
merge country year using "~/Documents/Projects/Global Inequality/ SWIID v3.1/SWIIDv3_1.dta", ///
	 keep(gini_net) _merge(_m1) nokeep
 
replace year=year-1 if _m1==1
sort country year
merge country year using "~/Documents/Projects/Global Inequality/ SWIID v3.1/SWIIDv3_1.dta", ///
	 keep(gini_net) _merge(_m1a) nokeep update
replace year=year+1 if _m1==1

sort country year
merge country year using "~/Documents/Projects/Data/GDPpc/pwt70_06032011version/PWT70_cysort.dta", ///
	keep(rgdpl) _merge(_m3) nokeep
replace rgdpl=rgdpl/1000

drop _m*
sort cc year
drop cc
egen cc=group(country)
sort cc year
saveold trend2.dta, replace

rsource using relig_chapter.R, rpath("/usr/bin/r") roptions(`"--vanilla"')

insheet using "trend3.csv", comma clear
capture drop v1
save trend3.dta, replace

use trend3.dta, clear
gen ch2=plfit
replace ch2=church2 if !mi(church2)
tsset cc year
xtreg ch2 l.gini_net l.rgdpl l.ch2, fe
xtreg gini_net l.gini_net l.rgdpl l.ch2, fe
xtreg rgdpl l.gini_net l.rgdpl l.ch2, fe


//Generate ten Monte Carlo simulations of the plfit series, replacing estimates with church2 if available						
set seed 3166							
forvalues a = 1/10 {
	gen e0 = rnormal()
	quietly tssmooth ma e00 = e0, weight (1 1 <2> 1 1)
	quietly sum e00
	quietly gen c`a'=plfit+e00*(1/r(sd))*plsefit
	replace c`a'=church2 if !mi(church2)
	drop e0 e00
}

save trend4.dta, replace


//Perform analysis using each of the ten simulations, saving the results
use trend4.dta, clear
local n_ivs = 4				
matrix coef = J(`n_ivs', 10, -99)
matrix se = J(`n_ivs', 10, -99)
matrix r_sq = J(3, 10, -99)
forvalues a = 1/10 {
	quietly xtreg c`a' l.gini_net l.rgdpl l.c`a', fe		
	matrix coef[1,`a'] = e(b)'
	matrix A = e(V)
	forvalues b = 1/`n_ivs' {
		matrix se[`b', `a'] = (A[`b',`b'])^.5
	}
	matrix r_sq[1, `a'] = e(r2_w)
	matrix r_sq[2, `a'] = e(r2_b)
	matrix r_sq[3, `a'] = e(r2_o)
}		

local cases = e(N)

svmat coef, names(coef)
svmat se, names(se)
svmat r_sq, names(r_sq)


//Display results across all simulations
egen coef_all = rowmean(coef1-coef10)

gen ss_all = 0
forvalues a = 1/10 {
	quietly replace ss_all = ss_all + (coef`a'-coef_all)^2
}

egen se_all = rowmean(se1-se10)
replace se_all = se_all + (((1+(1/10)) * ((1/9)*ss_all)))

gen t_all = coef_all/se_all
gen p_all = 2*normal(-abs(t_all))

egen r_sq_all = rowmean(r_sq1-r_sq10)

gen vars = " " in 1/`n_ivs'
local i = 0
foreach iv in l.gini_net l.rgdpl "l.ch2" "Constant" {
	local i = `i'+1
	replace vars = "`iv'" in `i'
}
mkmat coef_all se_all p_all if coef_all~=., matrix(res_all) rownames(vars)
di "Predicting church attendance"
matrix list res_all, format(%9.3f)
quietly sum r_sq_all in 1
local r2 = round(`r(mean)', .001)
di "R-sq within = `r2'"
quietly sum r_sq_all in 2
local r2 = round(`r(mean)', .001)
di "R-sq between = `r2'"
quietly sum r_sq_all in 3
local r2 = round(`r(mean)', .001)
di "R-sq overall = `r2'"
di "N = `cases'"

keep vars
svmat res_all, names(col)
drop if p_all==.
saveold ch2_res.dta, replace

use trend4.dta, clear

//Perform analysis using each of the ten simulations, saving the results
foreach v of varlist gini_net rgdpl {
	use trend4.dta, clear
	local n_ivs = 4				
	matrix coef = J(`n_ivs', 10, -99)
	matrix se = J(`n_ivs', 10, -99)
	matrix r_sq = J(3, 10, -99)
	forvalues a = 1/10 {
		quietly xtreg `v' l.gini_net l.rgdpl l.c`a', fe	/*to be replaced with your analysis*/	
		matrix coef[1,`a'] = e(b)'
		matrix A = e(V)
		forvalues b = 1/`n_ivs' {
			matrix se[`b', `a'] = (A[`b',`b'])^.5
		}
		matrix r_sq[1, `a'] = e(r2_w)
		matrix r_sq[2, `a'] = e(r2_b)
		matrix r_sq[3, `a'] = e(r2_o)
	}		
	
	local cases = e(N)
	
	svmat coef, names(coef)
	svmat se, names(se)
	svmat r_sq, names(r_sq)
	
	
	//Display results across all simulations
	egen coef_all = rowmean(coef1-coef10)
	
	gen ss_all = 0
	forvalues a = 1/10 {
		quietly replace ss_all = ss_all + (coef`a'-coef_all)^2
	}
	
	egen se_all = rowmean(se1-se10)
	replace se_all = se_all + (((1+(1/10)) * ((1/9)*ss_all)))
	
	gen t_all = coef_all/se_all
	gen p_all = 2*normal(-abs(t_all))
	
	egen r_sq_all = rowmean(r_sq1-r_sq10)
	
	gen vars = " " in 1/`n_ivs'
	local i = 0
	foreach iv in l.gini_net l.rgdpl "l.ch2" "Constant" {
		local i = `i'+1
		replace vars = "`iv'" in `i'
	}
	mkmat coef_all se_all p_all if coef_all~=., matrix(res_all) rownames(vars)
	di "Predicting `v'"  
	matrix list res_all, format(%9.3f)
	quietly sum r_sq_all in 1
	local r2 = round(`r(mean)', .001)
	di "R-sq within = `r2'"
	quietly sum r_sq_all in 2
	local r2 = round(`r(mean)', .001)
	di "R-sq between = `r2'"
	quietly sum r_sq_all in 3
	local r2 = round(`r(mean)', .001)
	di "R-sq overall = `r2'"
	di "N = `cases'"
	
	keep vars
	svmat res_all, names(col)
	drop if p_all==.
	saveold `v'_res.dta, replace
}

rsource using relig_chapter2.R, rpath("/usr/bin/r") roptions(`"--vanilla"')

shell open coef_graph.pdf

use trend4.dta, clear

// United States only
use trend4.dta, clear
local n_ivs = 4	
local n_dvs = `n_ivs' - 1 
local n_b = `n_ivs'*`n_dvs'/*VAR*/
matrix coef = J(`n_b', 10, -99)
matrix se = J(`n_b', 10, -99)
matrix r_sq = J(`n_dvs', 10, -99)
forvalues a = 1/10 {
	quietly var c`a' gini_net rgdpl if country=="United States", lags(1)		
	matrix coef[1,`a'] = e(b)'
	matrix A = e(V)
	forvalues b = 1/`n_b' {
		matrix se[`b', `a'] = (A[`b',`b'])^.5
	}
	matrix r_sq[1, `a'] = e(r2_1)
	matrix r_sq[2, `a'] = e(r2_2)
	matrix r_sq[3, `a'] = e(r2_3)
}		

local cases = e(N)

svmat coef, names(coef)
svmat se, names(se)
svmat r_sq, names(r_sq)

*display results
egen coef_all = rowmean(coef1-coef10)

gen ss_all = 0
forvalues a = 1/10 {
	quietly replace ss_all = ss_all + (coef`a'-coef_all)^2
}

egen se_all = rowmean(se1-se10)
replace se_all = se_all + (((1+(1/10)) * ((1/9)*ss_all)))

gen t_all = coef_all/se_all
gen p_all = 2*normal(-abs(t_all))

egen r_sq_all = rowmean(r_sq1-r_sq10)

// The rest needs work . . .
//gen vars = " " in 1/`n_ivs'
//local i = 0
//foreach iv in l.gini_net l.rgdpl "l.ch2" "Constant" {
//	local i = `i'+1
//	replace vars = "`iv'" in `i'
//}
//mkmat coef_all se_all p_all if coef_all~=., matrix(res_all) rownames(vars)
//di "Predicting church attendance"
//matrix list res_all, format(%9.3f)
//quietly sum r_sq_all in 1
//local r2 = round(`r(mean)', .001)
//di "R-sq within = `r2'"
//quietly sum r_sq_all in 2
//local r2 = round(`r(mean)', .001)
//di "R-sq between = `r2'"
//quietly sum r_sq_all in 3
//local r2 = round(`r(mean)', .001)
//di "R-sq overall = `r2'"
//di "N = `cases'"
//
//keep vars
//svmat res_all, names(col)
//drop if p_all==.
saveold us_res.dta, replace


//Single-country VARs
foreach cntry in "United States" Germany  {
	use trend4.dta, clear
	local n_ivs = 12				
	matrix coef = J(`n_ivs', 10, -99)
	matrix se = J(`n_ivs', 10, -99)
	matrix r_sq = J(3, 10, -99)
	forvalues a = 1/10 {
		quietly var gini_net rgdpl c`a' if country=="`cntry'", lags(1) /*to be replaced with your analysis*/	
		matrix coef[1,`a'] = e(b)'
		matrix A = e(V)
		forvalues b = 1/`n_ivs' {
			matrix se[`b', `a'] = (A[`b',`b'])^.5
		}
		matrix r_sq[1, `a'] = e(r2_1)
		matrix r_sq[2, `a'] = e(r2_2)
		matrix r_sq[3, `a'] = e(r2_3)
	}		
	
	local cases = e(N)
	
	svmat coef, names(coef)
	svmat se, names(se)
	svmat r_sq, names(r_sq)
	
	
	//Display results across all simulations
	egen coef_all = rowmean(coef1-coef10)
	
	gen ss_all = 0
	forvalues a = 1/10 {
		quietly replace ss_all = ss_all + (coef`a'-coef_all)^2
	}
	
	egen se_all = rowmean(se1-se10)
	replace se_all = se_all + (((1+(1/10)) * ((1/9)*ss_all)))
	
	gen t_all = coef_all/se_all
	gen p_all = 2*normal(-abs(t_all))
	
	egen r_sq_all = rowmean(r_sq1-r_sq10)
	
	gen vars = " " in 1/`n_ivs'
	local i = 0
	foreach iv in g_gini g_rgdpl g_ch2 "Constant" r_gini r_rgdpl r_ch2 "Constant" c_gini c_rgdpl c_ch2 "Constant" {
		local i = `i'+1
		replace vars = "`iv'" in `i'
	}
	mkmat coef_all se_all p_all if coef_all~=., matrix(res_all) rownames(vars)
	di "Predicting `cntry'"  
	matrix list res_all, format(%9.3f)
	quietly sum r_sq_all in 1
	local r2 = round(`r(mean)', .001)
	di "R-sq gini = `r2'"
	quietly sum r_sq_all in 2
	local r2 = round(`r(mean)', .001)
	di "R-sq rgdpl = `r2'"
	quietly sum r_sq_all in 3
	local r2 = round(`r(mean)', .001)
	di "R-sq ch2 = `r2'"
	di "N = `cases'"
	
	if "`cntry'"=="United States" local cntry="US"
		
	keep vars
	svmat res_all, names(col)
	drop if p_all==.
	saveold `cntry'_res.dta, replace /*"*/
}



set more on
log close
